#ifndef _toroidal_h_
#define _toroidal_h_

#define pi 3.1415926535
#include <stdio.h>
#include <gsl/gsl_spmatrix.h>

const int Nx = 192;
const int Ny = 256;
const int Nphi = 16;
const double Rmin = 3.069971;
const double Rmax = 5.069971;
const double Zmin = -1.495989;
const double Zmax = 1.495989;
const double dr = 0.010471;
const double dz = 0.011733;

const double rhoh_vs_a = 0.05; // Lamor radius / a
const double vh_vs_vA = 0.7; // fast particle velocity / Alfven velocity

const double psiaxis = -0.2443;
const double Raxis = 4.23;
struct vec
{
    double x, y ,z;
};
struct cylv
{
    double r,p,z;
};
struct pfvec
{
    vec E;
    vec B;
    vec gB;
    vec cb;
};
struct phsp // phase space (x,y,z,vparallel)
{
    vec X;
    double vp;
};

class vfield;
class vfiq;
class field
{
    public:
        field();
        field(double* dt);
        void load_fromfile(const char* str);
        double interp(cylv vp);
        vfield gradient();
        double* getdata() {return this->data; }
        void scale(double fact);
        void output(const char* str);
        void vtkoutput(const char* str);
    private:
        double data[Nx*Ny];

};

class vfield
{
    public:
        vfield();
        vfield(double* a, double* b, double*c );// (R, phi, Z)
        vfield(field* a, field* b, field* c);// (R, phi, Z)
        vec interp(cylv cord);
        void curl(vfield* vfp);
        void scale(double x);
        double* getfr(){return this->fr;}
        double* getfz(){return this->fz;}
        double* getfphi(){return this->fphi;}
        void vtkoutput(const char* str);
    private:
        double* fr;
        double* fz;
        double* fphi;
};

class particle
{
    public:
        particle() {}
        particle(double r,double z, double phi, double mu, double vpara);
        particle(cylv ccrd, double mu, double vpara);
        void push(vec E,vec B, vec gB,vec curb,double dt);
        cylv getcrd() {return this->crd;}
        double getvpara() {return this->vpara;}

        void setcrd(cylv crdn, double vpan) {this->crd = crdn; this->vpara=vpan;}
        void setmu(double mun) {this->mu = mun;}
        double Pphi();
        double energy(double B);
        double pmu(){return mu;}
        static double q;
        static double m;
    private:
        double vpara;
        cylv crd;
        double mu;
      
};

vec operator+(vec a, vec b);
vec operator*(double a, vec x);
double dot(vec a, vec b);
vec cross(vec a, vec b);
void print(vec a);

vec cyl2xyz(cylv a);
cylv xyz2cyl(vec a);
vec rotate(vec a, double phi); //vector field rotate

void partialy(field a, double* parx);
void partialx(field a, double* pary);

field fadd(field a, field b);
field fmul(field a, field b);

class world;

class ptcs
{
    public:
        ptcs(int N);
        ~ptcs();
        void init_hammersley(field B);
        void init_realistic(field B);
        void init_timeindependent(world wld);
        void output_realspace(const char* str);
        void output_velocityspace(const char* str);
        void output_poloidal(const char* str);
        double histo(double r);
        void scatter(double* pdata);
        double distribfunc(particle part);
        double pol_dstrb(cylv crd);

        int npt;
        particle* ppts;
        double* weights;
        double* periods;
        double* wdelta;
        field* pB;
        bool* inside;
        field* denfunc;

};
// in phase and quadrature
class vfiq;
class fiq;

class fiq
{
    private:
        double* f_inphase;
        double* f_quadrature;
    public:
        int n; // toroidal mode number
        fiq() {n = 0;}
        fiq(int n);
        fiq(double* a, double* b, int n);
        fiq(float* a, float* b, int n);
        fiq(field* a ,field* b, int n);
        double* get_fi();
        double* get_fq();
        double interp(cylv crd);
        void partialphi(fiq* pp);
        void gradient(vfiq* grad);
        void output(const char* str);
        void vtkoutput(const char* str);
        
};

class vfiq
{
    private:
        double* ifr;  // inphase component
        double* ifz;
        double* ifp;
        double* qfr;  //quadrature component
        double* qfz;
        double* qfp;
    public:
        int n; // toroidal mode number
        vfiq() {n = 0;}
        vfiq(int n);
        vfiq(vfield* a, vfield* b, int n);

        double* get_ifr(){return this->ifr;}
        double* get_ifz(){return this->ifz;}
        double* get_ifp(){return this->ifp;}
        double* get_qfr(){return this->qfr;}
        double* get_qfz(){return this->qfz;}
        double* get_qfp(){return this->qfp;}

        void curl(vfiq* vfiqp);
        void divergence(fiq* fiqp);
        vec interp(cylv crd);

};

bool checkboundary(cylv vp);
bool checkboundary_core(cylv vp);
class world
{
    public:
        field* Bsca;
        vfield* b;
        vfield* B0;
        vfield* gB;
        vfield* cb;
        field* vA;
        
        vfiq* E0;
        vfield* E1i;
        vfield* E1q;
        fiq* phies; // perturb electro potential
        fiq* denf; // kinetic compression density field

        gsl_spmatrix* partialr;
        gsl_spmatrix* partialz;
        gsl_spmatrix* bdotnabla;
        ptcs* fastptcs;
        field* fast_distrib;

        void init();
        void flb_gsl(fiq xiq, fiq yiq);
        void inertial_gsl(fiq xiq, fiq yiq);
        void shear(fiq phiiq, vfiq Bperpiq);
        void flow(fiq phiiq, vfiq viq);
        pfvec emfield(cylv crd);
        pfvec emfield_time(cylv crd, double tphase);
        phsp push_rk2(pfvec fs, particle* pt, double dt, double vpara);
        
        void partr_noiq(double* px, double* py);
        void partz_noiq(double* px, double* py);

        double wpi_single(particle pa, double omega);
        void wpi_multis_kc(double omega);
        void kcp_single(particle pa, double omega);

};

void perturb_hermitefunc(fiq* gen, int nh, int mp, double rho0, double shear, bool iq);
void trace_output(const char* str, cylv* pts, int n);
void perpgrad(fiq xiq, vfiq Eiq);
double innerproduct(fiq a, fiq b);
double inner_vf(vfiq F1, vfiq F2);

void eigen(double* data, int dim);
void eigen_gen(double* A, double* B, int dim, int eign, double* coefs);
void trace_vtkoutput(const char* str, cylv* pts, int n);
void linear_suppos(double* coefs, fiq* eigvecs, int n, fiq* eigm);

#endif